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We obtain the most general matrix criterion for stability and instability of multi- component soli- 
tary waves considering a system of N incoherently coupled nonlinear Schrodinger equations. Soliton 
stability is studied as a constrained variational problem which is reduced to finite-dimensional linear 
algebra. We prove that unstable (all real and positive) eigenvalues of the linear stability problem 
for multi-component solitary waves are connected with negative eigenvalues of the Hessian matrix, 
the latter is constructed for the energetic surface of N— component spatially localized stationary 
solutions. 



I. INTRODUCTION 



Recent discovery of self-focusing of partially coherent light and experimental observation of the so-called incoher- 
ent spatial solitons has called for a systematic analysis of the properties and stability of multi- component and 
multi-parameter solitary waves. Incoherent solitons are generated in noninstantaneous nonlinear media such as bi- 
ased photorefractive crystals. In this case, a self-consistent modal theory which is equivalent to the coherent 
density approach, describes the incoherent solitons with the help of a system of coupled nonlinear Schrodinger (NLS) 
equations (see also |§-pl). Similar models appear, in a different physical context, in the theory of soliton wavelength- 
division multiplexing M], multi-channel bit-parallel- wavelength optical fiber networks M, multi-species and spinor 
Bose-Einstein condensates [^), and other important applications H. In all such physical models solitary waves are 
multi- component, being described by localized solutions of the coupled nonlinear equations. In some very special cases, 
the coupled system allows for explicit analytical solutions (see, e.g., Ref. but, generally speaking, the nonlinear 
models with multi-component solitary waves are nonintegrable. Stability of solitary waves is therefore a crucial issue 
for any kind of their applications. 

The study of soliton stability has a long history. Stability of one-parameter solitary waves has been already well 
understood for both fundamental (single-hump and nodeless) solitons |lO|-|l^] and solitons with nodes and multiple 
humps HQ. The pioneering results of Vakhitov and Kolokolov |l0| found their rigorous justification in a general 
mathematical theory of Grillakis, Shatah, and Strauss Jig ]. Although the corresponding stability and instability 
theorems for the scalar NLS models extend formally to the case of multi-parametric solitons jlq] , all the cases 
analyzed so far correspond to solitary waves with effectively a single parameter. 

A recent progress in the study of soliton instabilities is associated with the application of a bifurcation theory valid for 
weakly unstable stationary localized waves. In this method, the corresponding unstable eigenvalue of the associated 
spectral problem is treated as a small parameter of multi-scale asymptotic expansions p6| . In the case of multi- 
parameter solitary waves, a simplified version of this method is usually reduced to a number of "magic determinants" 
constructed from the derivatives of the system invariants near a marginal stability line |l7| -|2C|] . However, such a 
bifurcation method has no rigorous proof, and it does not allow to predict the complete domains of the soliton 
stability and instability, since more general oscillatory instabilities may occur as well [ p^ , pl|p^ ] . 

In this paper, we present a complete theory for stability and instability of multi-parameter solitary waves considering 
a particular example of a system of N incoherently coupled NLS equations. Our results include the asymptotic 
bifurcation method with the determinant criterion as a simple near-threshold limiting case. They also expand the 
applicability boundaries of the previously known mathematical theorems fl5| to the case of multi-component and 
multi-parameter solitary waves. 

The system of incoherently coupled NLS equations has already been studied in many papers (see, e.g., Refs. j2^-|||] 
to cite a few). However, the study of stability of single- hump and multi- hump solitary waves was restricted again by 
a single-parameter case, when soliton components have a similar shape and their amplitudes are proportional to each 
other [ p3| . In this paper, we expand those results and present, for the first time to our knowledge, a complete matrix 
analysis of the constrained variational problem leading to the finite-dimensional linear algebra. Although some of our 
results depend on the properties which are specific to the model under consideration, we believe that both the method 
and the basic results can be generalized, under proper assumptions, to be applied to other types of nonlinear physical 
models that support multi-parameter solitary waves. 
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II. MODEL AND BASIC RESULTS 



We consider nonlinear interaction of N optical modes that describe either a multi-mode structure of a partially 
incoherent self-trapped beam or incoherent coupling between optical channels with different wavelengths in a fiber. 
Then, the amplitude envelopes of partial modes satisfy the following system of incoherently coupled NLS equations, 

.dipn + ^ n Y^ n _j_ I \ " j nm \tp m \' 2 ] tj) n = 0, (1) 



dz 



\m=l ) 



where stands for Laplacian in the D— dimensional space x = [x\, ...,xd), and all the coefficients d n are assumed 
to be positive. When one of the variables of the vector x stands for time, Eqs. (|l|) describe the spatio-temporal 
dynamics of self-focused and self-modulated light in the form of the so-called light bullets. 

Provided the symmetry conditions j nm = Jmn are satisfied, the system (Q) conserves the Hamiltonian, 

,oo / n n n \ 

H= dx K]4 |V X ^„| 2 -"^^ Tnml^PlV'ml 2 , 

J -° c \n=l „=l m =l / 

the individual mode powers, Q n = 5 / |^n| 2 rfx, and the total field momentum. Localized solutions of Eqs. (^) for the 
fundamental solitary waves are defined as ip n = $ Tl (x)e l ' 3 " z , where $ n (x) are real functions with no nodes, and /3„ 
are positive propagation constants. The soliton solutions are stationary points of the Lyapunov functional, 

JV 

A[^] = H[^} + J2^Q n m (2) 

n=l 

i.e. the first variation of A[i/>] vanishes at ip = $(x). The second variation of A[i/>] defines the stability properties: 
negative directions of the second variation corresponds to unstable eigenvalues in the soliton stability problem (see, 
e.g., Ref. for a review of the basic results). 

The stability problem is defined by minimizing the second variation of the Lyapunov functional A [■«/>], 

/oo 
dx[(u|Liu) + (w|L w>], (3) 
-00 

where u(x) and w(x) are perturbations of the multi-component solitary wave taken in the form ip = 3>(x) + 
[u + iw] (x)e Az , and the scalar product is defined as (f|g) = Y^n=i fn9n- The matrix Sturm-Liouville operator 
Lq has a diagonal form with the elements 



N 

- m ) 

771 — 1 



(L )nn = -dn^l +Pn~Y^ Inm^l 

and the matrix operator Li has the elements 



N 



m — 1 



at the diagonal, and {L\) nm — — 27„ m $„$ m , off the diagonal. Operators Lq and L\ determine the linear eigenvalue 
problem for stability of multi-component solitary waves, 

Liu = —Aw, Low = Au. (4) 

Both the linear problem (Q) and minimization problem (||) should satisfy a set of N constraints, 



F n = I dx($„e„|u) = 0, (5) 
■respc 

of a perturbation described by a vector (u,w). 



where e„ is the unit vector, which correspond to the conservation of the individual powers Q n under the action 
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First of all, we recall the main result of Refs. [|l 0| |l 2|] that the one-parameter solitary waves with no nodes (N = 1) 
are stable in the framework of the constrained variational problem (|J)-(||) provided the energetic surface A s (/3) = A[<&] 
is concave up, i.e. 

d 2 A s dQi f . 

Under this condition, the linear eigenvalue problem (^) has no unstable eigenvalues, i.e. those with a positive real 
part A. Otherwise, the second variation (||) constrained by the set (|J) has a single negative direction that corresponds 
to a single positive eigenvalue A in the linear eigenvalue problem (||) jl^Jl^]. The stability criterion for scalar (or 
one-component) NLS solitons holds when the self-adjoint operator Li has a single negative eigenvalue, i.e. when the 
second variation (|^), without the constraint (|J) imposed, has a single negative direction. If the latter condition is 
not satisfied, as it happens for solitary waves with nodes, the fundamental criterion for the soliton instability can be 
extended only for a special case p^O], while more generic mechanisms of oscillatory instabilities, associated with 
complex eigenvalues of the linear eigenvalue problem, may appear beyond the prediction of the fundamental criterion 

Here we extend the soliton stability analysis to the case of multi-component solitary waves described by a system 
of incoherently coupled NLS equations ([!]) . We assume that the number of negative directions (eigenvalues) of the 
second variation S 2 A is fixed, and we denote it as n(A). The unstable eigenvalues A of the linear problem (|]) are 
connected with some negative eigenvalues of the matrix U defined by the elements 

d 2 A s _ dQ n _ dQ m , s 

nm ~ dp n 8(3 m - df3 m d/3 n - U 

The matrix U is the Hessian matrix of the energetic surface A s ((3) . We denote the number of positive eigenvalues 
of the matrix U as p(U), and the number of its negative eigenvalues as n(U), so that p(U) + n{U) < N, since some 
eigenvalues may be zeros in a degenerate (bifurcation) case. As is shown below, both p(U) and n(U) satisfy some 
additional constraints, 

p(U) < min{7V, n(A)}, n(U) > max{0, N - n(A)}. (8) 

Within these notations, we formulate (and prove below) the following fundamental results on stability and instability 
of multi-component solitary waves of the coupled NLS equations (^) : 

(i) the linear problem (^) may have at most n(A) unstable eigenvalues A, all real and positive; 

(ii) a multi-component soliton is linearly unstable provided p(U) < n(A); then the linear problem (^) has n(A)— p(U) 
real (positive or zero-becoming-positive) eigenvalues A; 

(iii) a multi-component soliton is linearly stable provided p(U) — n(A)(< N); in the case n(A) = N this criterion 
implies that the energetic surface A s (/3) is concave up in the /3-space; 

(iv) a single eigenvalue A crosses a marginal stability curve when the matrix U possesses a zero-becoming-negative 
eigenvalue; the normal form for the instability-induced dynamics of multi-component solitary waves resembles the 
equation of motion for an effective classical particle subjected to a iV-dimensional potential field, 

1 N N A A 

n— 1 m— 1 

where M nm are the elements of the positive-definite "mass matrix" [see Eq. ( |2S| ) below] , v is the vector describing a 
perturbation to the soliton parameters (3, and W((3, v) is an effective potential energy defined as 

N 

W((3, v) = H s ((3 + v) - H s {f3) + {fin + v n ) [Q sn (P + u)- Q m {fi)] ■ (10) 

n=l 

These results should be compared with the results following from the stability and instability theorems earlier 
formulated by Grillakis, Shatah, and Strauss |H| ]. The stability result (iii), i.e. the condition p(U) = n(A), is identical 
to that of the stability theorem M|, but the instability results (i)-(ii) are more general and explicit. In particular, 
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the theorem of Grillakis et al. guarantees the soliton instability provided the difference n(A) — p{U) is odd. 
However, our results predict that the soliton instability always occurs for n(A) — p(U) > 0, being associated with 
exactly n(A) — p(U) non-negative real eigenvalues A of the linear eigenvalue problem (^). Moreover, according to 
our result (iv), each new unstable eigenvalue A appears via a bifurcation at the marginal stability curve where the 
determinant of the matrix U vanishes, i.e. it is connected with a zero-becoming-negative eigenvalue of the Hessian 
matrix U. If n(A) > N, unstable eigenvalues A originated from the negative eigenvalues of the Hessian matrix U 
co-exist with n(A) — TV unstable eigenvalues of the linear problem (|J), i.e. a solitary wave is unconditionally unstable 
when n(A) > N. 



III. A PROOF OF THE BASIC RESULTS 



Here we develop the analysis of the problem (|3|)-(||), in order to prove the results (i)-(iv) formulated above. The 
Sturm-Liouville operators {L ) nn are all non- negative since the fundamental (node- less) solution <&(x) for a solitary 
wave reaches a bottom of the spectrum at zero: (Lo)nn$n = 0. As a result, 



(min) 5 2 A = \ dx (ujlgu) = / dx (u fc |u fc ). (11) 



Here (//,Ufc) are eigenvalues and eigenfunctions of the auxiliary linear problem, 



JY 



LiUfc = ^iUfe - ^2 ^m$m(x)e m . (12) 
m— 1 

The linear problem (|l2] ) is constrained by the set (||) and the parameters V\, V2,...,vn have the meaning of the Lagrange 
multipliers. 

Let us suppose that the Sturm-Liouville matrix operator Li has n(A) negative eigenvalues /i = 
{M-n(A)>M-n(A)+i>—)M-l} corresponding to the eigenfunctions u = {i/>_ n(A) (x), ip^ n{A)+1 (x.), ■0_ 1 (x)}; a sin- 
gle zero eigenvalue with a one-node eigenfunction u = d^/dx and the rest of the spectrum is positive and con- 
tains N branches of the continuous spectrum, for /i > /?2, /3jv}, and some isolated positive eigenvalues, for 
i" = {Mi,M2, ■■■,A 1 p}- 

The mathematical problem can then be reformulated in the following way. The linear operator Li has n(A) negative 
eigenvalues that generate negative directions of the second variation S 2 A. However, the corresponding eigenfunctions 
do not satisfy generally the constraints By introducing the Lagrangian multipliers in (|ll|) and ([l2]), we satisfy 
a constrained minimization problem (Q) and ^) but, due to this procedure, the number of negative eigenvalues may 
be reduced. We will show how to connect the total number of negative eigenvalues of the constrained problem 
(HI), and (H) with the negative eigenvalues of the Hessian matrix ffl). But, as a prerequisite, we prove two additional 
results for the spectrum of the problem (||): 

(i) the spectrum of A 2 is real, i.e. oscillatory instabilities are prohibited; 

(ii) each negative direction (/i, Ufc) of the problem (|l2| ) generates an unstable (positive) eigenvalue A of the problem 
(§• 

To prove the statement (i), we notice that the matrix operator Lo can be factorized as Lo = Ylct=i M^M^, where 
Md has a diagonal form with the following matrix elements 



(M d ) nn = 



provided the soliton solutions < i ) n (x) have no nodes in a finite domain. Using this factorization, the linear problem 
(f|) can be rewritten for the function u = J2d=i v d as follows, 



D 



MdLxM+Vd, - -A 2 v d . 

d'=l 

Since the matrix operator with the elements M^LiM^, is Hcrmitian, its eigenvalues (—A 2 ) are all real. 
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To prove the statement (ii), we suppose that we have constructed a negative direction (/z, U&) of the problem 
l|)-(|l2l) subject to the constraints (||). Then, the linear problem (Q) has an unstable eigenvalue A defined as, 



2 _ (u fc |L Liu fc ) (u fe |L u fe ) 

A — j r — — fl — j r— . [L6j 

Since the linear operator Lo is positive definite for any uj- 7^ (0, <J>), we have A 2 > for any /i < 0. 

Our next goal is to construct solutions to the auxiliary problem ([l2]). Since the matrix operator Li is Hermitian, 
it has a complete spectrum in a Hilbert space that is suitable for expanding the vector function u&(x). We present 
such a spectral decomposition in the form, 



<W = 2^™ 2^ u _ u — */v( x ) + 

m=l Vr<0 ^ ^ r ft r >0 



u,.lx) > ; r,„ [ > V',(x) 4- r ^lg^ ^(x)) , (W) 

M,->0 A* Mr / 



where the sum <0 contains n(A) terms from the negative spectrum, while the sum >0 includes schematically 
both the discrete and continuous positive spectra of Lj. The contribution from the neutral eigenfunction u = d^/dx 
vanishes due to the symmetry properties. General solution (|l^) has to be constrained by the conditions (|J). This 
system reduces to the linear algebra for the Lagrange multipliers, A.(fi)v = 0, where the matrix A(/x) has a symmetric 
form with the elements: 

A„ m (M) = V + V ■ (15) 

i< — u r ^— a — a r 

ti r <o ^ ^ r ti r >o ^ ^ r 

The linear system A.(fi)u = 71/ has generally TV real eigenvalues 71 (/x), 72(^)1 ••■7 7jv(m)- These eigenvalues are 
continuous functions of for < 0, except for n(A) resonant planes at fi = {/i-n(A)> M— n(A)+i) •••! A*-i}- At these 
planes, the matrix A(/u) has poles and the eigenvalues y(fi) may have singularities. Below, we prove the following 
three properties of the eigenvalues y(fi): 

(i) all eigenvalues 7(/k) are negative for /1 < A*_ n (A)(< 0); 

(ii) each eigenvalue j(fi) is a decreasing function of fx for /i < 0, except for n(A) resonant planes at // = 

U*-n(A)> M-n(A)+lj — 

(iii) at least (TV — 1) eigenvalues 7(/i) are continuous at any of the resonant planes fi = \i r < 0, while the minimal 
eigenvalue, say 71 (/i), may have an infinite discontinuity jumping from negative infinity, at jtt — * fi r — 0, to positive 
infinity, at fi — ► /!,. + 0. 

To show the property (i), we consider the asymptotic limit of A(fi) as — » —00. In this limit, the eigenvalues j(fi) 
can be expressed from the algebra of quadratic forms, 



X 1 ' Vr<0 U r >0 / 



7(M) = — 7— ( > . '' - > , ) • (.Ui) 

where 



6 r = 



jv 2 
^ i/ n (t/v|$„e„ 



> 0. (17) 



Since all b r may not vanish simultaneously for u ^ 0, the eigenvalues j(fi) are negative definite in Eq. ( |l6| ) so that 
7(/i) — ► —0 as /i — > —00. 

To show the property (ii), we take the derivative of the system A.u = 'j(fi)i f and use the algebra of quadratic forms. 
The derivative of j(fi) is then defined for [i < excluding the resonant planes at // = {/i_ n (A), /i- n (A)+ij 1} as 



e^M = 1 , , dA(/x) 1 

d/i tfyi (fl - fl r ) 2 ' ^ Q (fl - fl r ) 2 
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where b r are defined by the same relation (^). Since the derivative of j(p) is hegative definite in Eq. (|T^), all 
eigenvalues 7(/z) are decreasing functions of \i whenever dhf(p)/dp exists. 

To show the property (iii), we consider the behavior of the eigenvalues 7(ju) at the resonant plane fi — fj, r < 0. In 
this limit, the matrix elements A nm (fi) have the following asymptotic form, 

. ( , ($„e n |Vv) (tp r \<S> m e m ) 



(ji - fl r ) 

Therefore, the matrix A(/z) has (N — 1) zero eigenvalues "f(p) and a single non-zero eigenvalue 71 (/x) with the 
asymptotic value, 

1 * 

7i(M)-7— pr^l^enWI 3 . (19) 
Mr J n=1 

If the sum in Eq. ( |l9| ) does not vanish, the eigenvalue 71 (/i) has an infinite discontinuity described in (iii) and, 
according to the property (ii), it is the minimal eigenvalue. Other (N — 1) eigenvalues are in fact non-zero in the 
limit /i — > fi r . Since the matrix A(/x) is a meromorphic function of fi as /1 < 0, the eigenvalue (fi — fi r )j(fi) is of order 
of 0(fi — fi r ) for (N — 1) non-singular eigenvalues. Therefore, the values of 7(/x) are generally non-zero in the limit 
/i -> /i r . 

Thus, we have a clear picture how the eigenvalues 7(/i) behave as functions of \x [see Figs l(a,b)]. Starting with 
small negative values as /i — + —00, all eigenvalues decrease as fi grows towards the n(A) resonant planes. At each of 
those planes, (N — 1) eigenvalues remain continuously decreasing, while one (minimal) eigenvalue jumps to a positive 
domain unless the condition 

N 

£|($„e„h/v)| 2 =0 (20) 

71=1 

is satisfied (this condition will be discussed below). Assuming that the condition ( pp| ) is not met, we come to the 
conclusion that a root of 7(/i) may occur only after a jump of j(fi) at a resonant plane fi = fi r to a large positive 
value, and further decrease of 7M as fi (> fi r ) grows. The root of j(fJ>), if exists for fi < 0, produces a legitimate 
solution Ufc(x) of the problem ( fl2|) under the constraints (||). This solution (/i, Ufc) would then be associated with an 
unstable eigenvalue A, according to the connection formula (p"3|). Thus, our main task is to control the behavior of 
positive 7(/x) in between the plane \i = and the resonant planes /1 = {^-n(A)j A t -n(A)+i: •••jM-l}- 
At the plane /i = 0, the problem ( |l2] ) has a simple solution for Ufe(x), 

u M=0 (x) = 5> n -^. (21) 

n— 1 

Substituting Eq. ( |2l] ) into the constraints (||), we find that A(0) = U, where the matrix U is the Hessian of the 
energetic surface A s (f3) with the elements U nm defined by Eq. (ffl). We can now use this construction and prove the 
main results (i)-(iii). In the analysis below we assume that the condition (2w is never met and the root of 7(/i) at 
(i < is associated with the unstable eigenvalue A of the stability problem (3). 

The roots of j(fJ-) may appear only to the right of any of the n(A) resonant plane. There are totally n(A) jumps of 
7(/«) to positive values at /i < and, therefore, no more than n(A) roots of 7(/x) may exist for /1 < 0. 

If n(A) = N, the positive eigenvalues 7(/i) remain continuous after passing the corresponding resonant plane at 
fi = fi r (< 0). Therefore, the sign of these eigenvalues is controlled by the eigenvalues of the Hessian matrix U at 
/i = 0. If p(U) = N = n(A), all positive eigenvalues 'y(fi) remain positive for \i r < /1 < and no roots of 7(/x) exist 
for /i < [see Fig. 1(a)]. If p(U) < N = n(A), there exist N — p(U) negative or zero-becoming- negative eigenvalues 
of U that correspond to N — p{U) roots of 7(//) for p, < 0. 

If n(A) < N, then N — n(A) eigenvalues 7(/i) do not have jumps at the corresponding resonant planes /1 = p, r . 
They continue to be negative and match at fi = with the N — n(A) negative eigenvalues of U. From this, we come 
to the conclusion that p(U) and n(U) satisfy the constraints i.e. n(U) > N — n(A) or, equivalently, p(U) < n(A). 
Furthermore, other n(A) (< N) eigenvalues 7(/x) may have roots for /j, < which are completely controlled by the 
remaining n(A) eigenvalues of U according to the same criterion as in the case p(A) = N. For instance, if p(U) < n(A), 
then n(A) —p(U) negative or zero-becoming-negative eigenvalues of the matrix U correspond to n(A) —p(U) roots of 
7(/i) at fi < 0. 
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If n(A) > N, then n(A) — N eigenvalues j(fi) jump twice in the domain fi < leading to at least n(A) — N 
unconditional roots for fi < [see Fig. 1(b)]. After the jumps, the N eigenvalues 'y(fi) match the N eigenvalues of the 
matrix U and may have additional roots of 7(/x) if p(U) < N. Total number of roots of jifi) at < is then defined 
as (n(A) — N) + (N — p{U)) = n(A) - p(U). 

The analysis above is valid for the non-degenerate case when the condition ( p0| ) is never satisfied. However, the 
stability and instability results (i)-(iii) are not affected even if the condition ( ^o|) is satisfied for a particular resonant 
plane fi = fi r (< 0). In this case, the eigenfunction u&(x) of the operator Li satisfies all the constraints dq ) identically 
and, therefore, the eigenvalue fi = fi r is associated with an unstable eigenvalue A, according to Eq. Jig). Although 
the eigenvalue 71 (/x) has no jump at fi = fi r [see Eq. (19)] and is continuous, it is still controlled by the negative 
eigenvalues of U at fi = 0. Indeed, in this case, the minimal eigenvalue Jx(fj) at fi < fi r remains negative for fi > \i r 
and matches with a negative eigenvalue of the matrix U (if no other jumps occur in the domain fi < 0). This additional 
negative eigenvalue fi still predicts the instability, according to the result (iii). 

Finally, we prove the result (iv) for the instability bifurcation of multi-component solitary waves. Provided the 
number n(A) is fixed, the instability bifurcation may only occur when A(0) = U has a zero eigenvalue for a certain 
eigenvector v = . Let us define U = Uthr at the marginal stability curve so that the determinant of XJthr vanishes. 
The instability bifurcations of multi-component solitons were considered in Refs. Jl7| , p0|]_ but the results do not agree 
with each other. Here, we recover the results of Ref. Jl?] ] and derive the normal forni(p|) by an elegant reduction of 
general algebraic expressions. 

Assuming fi = for i> — i/W so that \5 thr u^ = 0, we find the asymptotic solution of Eq. (Q) in the form (|2l| ) and 



w, = „(x)=Af;^L -« 



(22) 



In this limit, the second variation <5 2 A of the Lyapunov functional can be found from Eqs. (||), (pi]), and ( |22| ) as 
follows 



£ 2 A = D 1 X 2 



(23) 



where 



f/x 



N 



' N 

n=l 



/ X dx'$ m (x') 
Jo 



(24) 



The integral converges under the condition that is a solution of the equation Uthr^^- 1 = 0. On the other hand, 
the perturbation (21) shifts the soliton parameter (3 according to the expression, 3>(x, /3) + u M= o( x ) - * /3 + u^ k '). 



As a result, the second variation can be closed as 



5 2 A 



A- A, + 



(fe)- 



-D Q , 



(25) 



where 



D a = (i/ fc )|Ui/W) 



The parameter A in Eq. (|25|) is chosen from the condition that the first variation of A s (/3 + u 

i/(fc). This gives the connection formula: A = A sf = H s ((3) + ^n=i(^" ' " 
(po]), we recover the result of the bifurcation theory, 



(26) 

' fc - ) ) vanishes for arbitrary 
)Q m (P). Equating Eq. (||) and Eq. 



A 



Do 



(27) 



Since Di > [see Eq. (|I])], the positive values of A 2 occur when the determinant of the matrix U is small and negative 
(i.e. the matrix U has a zero-becoming-negative eigenvalue when the soliton parameter f3 crosses the marginal stability 
curve). The explicit formulas of the soliton bifurcation theory provide an alternative and more compact form for the 
determinants Do and D\ compared to those obtained in Ref. |20|. 

The normal form (||) follows from Eqs. (|2^) and (^5|) when A = A st + E, and the perturbation vector is 
replaced by a slowly varying vector v = i/(zj(scc |l7j for details). Then, the surface A s (/3 + v) is extended beyond 
the second variation limit, and the linear approximation is converted into the slope: \i> (z) = du(z)/dz. The mass 
constants M nm follow from Eq. (p4|) in the explicit form: 
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- iy±m or w 2 ^) or ■ 

The normal form @ resembles the conserved sum of the kinetic energy and potential energy W(f3, v) of a particle 
moving in a iV-dimensional space. We notice that the kinetic energy with the "mass" matrix (^) is positive definite 
and the unperturbed multicomponent solitary wave (i.e. that with v = 0) is a stationary point of W(f3, u) for any 
v. Thus, the stability of multi-component solitons resembles the stability of a particle located at an equilibrium 
point of the iV-dimensional field p7| . Under the condition that n(A) = N, the particle is stable if in the /3-space 
the potential energy surface W{(3) is concave up, and it is unstable if the potential energy surface is saddle-type or 
concave down. If n(A) < N, the potential energy surface W(f3) always has some N — n(A) negative directions that 
do not affect the stability properties of the particle. However, the remaining n(A) (< N) directions of the potential 
energy surface define the stability of the particle with the same criterion as above. Finally, for the case n(A) > N, 
the soliton stability properties defined by the type of the potential energy surface W(f3) are not conclusive since the 
corresponding unstable eigenvalues coexist with additional n(A) — N unconditionally unstable eigenvalues. 



IV. EXAMPLE: TWO COUPLED NLS EQUATIONS 



In order to demonstrate how our general theory can be applied to a particular physical problem and also to compare 
the stability and instability results (ii)-(iii) with some earlier known exam ples , we consider here an important case of 
two incoherently coupled NLS equations in (1+1) dimension (see, e.g., |23 25|): 



.dtjj! 
oz 

i— 

dz 



dx 2 
dx 2 



(l^il 2 



+ (|^2| 2 +7l^l| 2 )V>2 = 0, 



(29) 



where 7 is a coupling parameter. System (|29|) is a two-component reduction of the general N— component system (|1|) 
for d\ = 0I2 = 1, 711 = 722 = 1) arid 712 = 721 = 7- An explicit soliton solution can be easily found for fa = fa = f3 
and 7 > — 1 in the form, 




sech (^ffix^ 



(30) 



This solution describes a two-component solitary wave with the components of equal amplitudes. It corresponds to a 
straight line fa = /3% in the parameter plane (fix, fa) of a general two-parameter family of solitary waves of the model 
(p9|). When — 1 < 7 < 0, such two-parameter solitons may exist everywhere in the plane (fa, fa), while for 7 > 0, the 
soliton existence domain is restricted by two bifurcation curves, fa = uj±(-y)fa, where 



w±(7) 



/ VTT8t"- 1 



±2 



(31) 



Approximate analytical expressions can also be obtained in the vicinity of the bifurcation curves @, when one 
of the components of a composite solitary wave becomes small, while the other one is described by the scalar NLS 
equation. Such a case, when one of the component creates an effective waveguide that guides the other component, is 
known to describe the so-called shepherding effect where the large-amplitude component plays a role of a shepherding 
pulse 0| . The composite soliton, that describes a shepherding pulse tpi guiding a small component pulse tp2, can be 
found in the form (see also Ref. p5[), 



$1 = R (x) + e 2 R 2 (x) + 0(e 4 ), $ 2 = eSx(x) + 0(e 3 ). 
It exists in the vicinity of the bifurcation curve, 

fa = w+( 7 )/3i + e 2 ^ 2 +(7)/3i + 0(e 4 ), 



(32) 



(33) 



and the main terms of the asymptotic series (32), (|33| ) are defined as 

R = V2/3isech (y/jhx) , #1 = VT^sech^ {y/fax) 
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and 



W2+ 



I- oc dx(Sf + gTgoggg) 



The second-order correction Rz{x) is a solution of the differential equation, 



-<9 2 + /3i - 6/3i sech* 



B.2 = 7fio5J. 



From the domain of existence of the two-component soliton, it follows that ^+(7) > 0, for < 7 < 1, and 
u, 2+{~f) < 0, for 7 > 1. At 7 = 1 (the so-called integrable Manakov case), a family of two-parameter composite 
solitons becomes degenerated: it exists on the line Pi = P 2 but, generally, it is different from the one-parameter 
solution (^0|). The coupled solitons are known to be stable for the integrable case 7 = 1. Here we apply the stability 
theory developed above and prove that the (l+l)-dimensional two-parameter solitons, including the solitons of equal 
amplitudes (pO|), are stable for 7 > 0, and unstable for 7 < 0. 

First, we evaluate the indices p(U) and n(A) for the explicit solution ([50j). As follows from Eqs. j29| ) and (|30|), the 
Hessian matrix U with the elements (0) can be found in the form, 



dQi dQ 2 



dpi 



dp 2 V^(l + 7) 



and 



dQ 1 = dQ 1 



V^(l + 7)' 



As follows from these results, the Hessian matrix has p(U) = 2 positive eigenvalues, for —1 < 7 < 1, and p(U) = 1 
positive eigenvalue, for 7 > 1. On the other hand, the linear matrix operator Li given below Eq. (§ can be 
diagonalized for linear combinations of the eigenfunctions, V\ = U\ + u 2 and v 2 — U\ — u 2l such that 



(3-6(3 sech 2 (jpx^j 



-d 2 x + 0-20^-^1 sech 2 (^ 3 



vi = fjtvi, 
v 2 = p,v 2 . 



(34) 



Both the operators in Eqs. (|34|) are the linear Schrodinger operators with solvable sech-type potentials, and the 
corresponding eigenvalue spectra are well studied. The first operator always has a single negative eigenvalue for 
fi = —3/3, whereas the second operator has no negative eigenvalues, for 7 > 1, has a single negative eigenvalue, for 
< 7 < 1, and it has two negative eigenvalues, for — 1 < 7 < 0. Thus, in total there exist n(A) = 3 negative 
eigenvalues, for — 1 < 7 < 0, n(A) = 2 negative eigenvalues, for < 7 < 1, and n(A) = 1 negative eigenvalue, for 
7 > 1. 

Applying the stability and instability results (ii)-(iii) obtained and discussed in Sees. II and III, we come to the 
conclusion that the soliton solution (^0|) with equal amplitudes is linearly stable for 7 > 0, since in this domain 
p(U) = n(A) = {1, 2}, and it is linearly unstable for — 1 < 7 < 0, since in this domain p(U) — 2 < n(A) = 3. 

The soliton stability in the model ( p9| ) for 7 > was also studied by Berge |^3j who considered the case of degenerate 
one-parametric solitary waves (30). Here we have extended those results to a general case: the same stability and 
instability results are valid for the two-parameter family of solitons provided that the indices p(U) and n(A) remain 
unchanged for the values {(3\,(3 2 ) that do not belong to the line (3\ — (3 2 . Indeed, applying a perturbation theory for 
small 7 (see Q for details), one can show that n(A) = 3, for — 1 <C 7 < 0, and n(A) = 2, for < 7 1. 

Finally, we consider the other limiting case that describes the shepherding effect [see Eqs. (^2|) and (|33|)]. In this 
limit, the elements (w of the Hessian matrix U can also be calculated in explicit analytical form, 



ggi 

d(3i 



1 



+ 0(e 2 ), 



dQi _ dQ 2 _ j2_ 

d(3 2 dpi uj 2+ 



+ 0(e 2 ), 



dQ 2 

d(3 2 



Sl 

u 2+ 



where 



si 



Sfdx and r 2 = 



RoR 2 dx. 



Since s\ > for any 7, while uj 2+ (j) > for < 7 < 1 and u 2+ {^) < for 7 > 1, the Hessian matrix U calculated 
for the shepherding soliton ( |32| ) has p(U) = 2 positive eigenvalues, for < 7 < 1, and p(U) — 1 positive eigenvalue, 
for 7 > 1. 
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On the other hand, the linear matrix operator Li cannot be diagonalized for the shepherding soliton (|32|) unless 
e = 0. In the latter (decoupled) case, it has a single negative eigenvalue at n = —3/3i and a double degenerate zero 
eigenvalue. When e / 0, the zero eigenvalue shifts to become fx = — 2w 2 +(7)/3ie 2 + 0(e 4 ). Therefore, the matrix 
operator Li for the shepherding soliton ( j32|) has n(A) = 2 negative eigenvalues, for < 7 < 1, and n(A) = 1 
negative eigenvalue, for 7 > 1. Thus, we come to the conclusion that the shepherding soliton is stable for 7 > since 
p(C0=n(A) = {l,2}. 



V. CONCLUSION 



We have developed a rigorous stability analysis of multi-component solitary waves considering a system of incoher- 
ently coupled NLS equations (0) as a particular but important physical example. The method and, correspondingly, 
both stability and instability results can be extended to other types of solitary waves, such as multi-component spatial 
solitons (e.g., incoherent solitons) in non-Kerr (e.g., saturable) media, parametric solitary waves in quadratic (or x*- 2 -*) 
optical media, etc. In all such cases, our stability and instability results (i)-(iv) can be readily generalized with a 
rigorous proof of some of the previously known results of the asymptotic multi-scale expansion theory. However, 
additional analysis is required in each of those cases in order to clarify the conditions when these results completely 
define the stability properties of multi-component solitary waves. In the cases beyond these conditions, oscillatory 
instabilities may occur, and appropriate studies should rely solely on the numerical analysis of the corresponding 
eigenvalue problems and their linear spectra. 



ACKNOWLEDGMENTS 



The authors appreciate collaboration with A. A. Sukhorukov and E.A. Ostrovskaya. Dmitry Pelinovsky thanks 
C.K.R.T. Jones, K. Promislow, M. Weinstein, J. Yang, and A. Yew for stimulating discussions at different stages of 
the work. The work is supported by the Performance and Planning Fund of the Institute of Advanced Studies and 
by the Australian Photonics Cooperative Research Centre. 



t Permanent address: Department of Mathematics, McMaster University, Hamilton, Ontario L8S 4K1, Canada. 
[1] M. Mitchell, Z. Chen, M. Shih, and M. Segev, Phys. Rev. Lett. 77, 490 (1996); M. Mitchell and M. Segev, Nature (London) 
387, 880 (1997). 

M. Mitchell, M. Segev, T.H. Coskun, and D.N. Christodoulides, Phys. Rev. Lett. 79, 4990 (1997). 

See, e.g., D.N. Christodoulides, T.H. Coskun, M. Mitchell, and M. Segev, Phys. Rev. Lett. 80, 2310 (1998); D.N. 
Christodoulides, T.H. Coskun, M. Mitchell, Z. Chen, and M. Segev, Phys. Rev. Lett. 80, 5113 (1998). 
Y. Nogami and C.S. Warke, Phys. Lett. A 59, 251 (1976); A. A. Sukhorukov and N.N. Akhmediev, Phys. Rev. Lett. 83, 
4736 (1999); N.N. Akhmediev, W. Krolikowski, and A.W. Snyder, Phys. Rev. Lett. 81, 4632 (1998). 
O. Bang, D. Edmindson, and W. Krolikowski, Phys. Rev. Lett. 83, 5479 (1999). 
See, e.g., S. Chakravarty, M.J. Ablowitz, JR. Sauer, and R.B. Jenkins, Opt. Lett. 20, 136 (1995). 

C. Yeh and L. B ergman, Phys. Rev. E 57, 2398 (1998); 60, 2306 (1999); Yu.S. Kivshar and E.A. Ostrovskaya, arXiv: 



patt-sol/9912008| , submitted to Opt. Lett. (2000). 

T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998); MR. Matthews et al, Phys. Rev. Lett. 83, 3358 (1999); S.-K. Yip, Phys. Rev. 
Lett. 83, 4677 (1999). 

A. Hasegawa, Opt. Lett. 5, 416 (1980); FT. Hioe, Phys. Rev. Lett. 82, 1152 (1999); B. Tan and J.P. Boyd, Chaos, Solitons 
and Fractals 11, 1113 (2000). 

N.G. Vakhitov and A. A. Kolokolov, Radiophys. Quantum Electron. 16, 783 (1973); A. A. Kolokolov, Radiophys. Quantum 
Electron. 17, 1016 (1976). 

E.A. Kuznetsov, A.M. Rubenchik, and V.E. Zakharov, Phys. Rep. 142, 103 (1986); V.G. Makhankov, Yu.P. Rybakov, and 
V.I. Sanyuk, Physics-Uspekhi 62, 113 (1994). 

J. Shatah and W. Strauss, Comm. Math. Phys. 100, 173 (1985); M.I. Weinstein, Comm. Pure Appl. Math. 39, 5 (1986). 

C. K.R.T. Jones, J. Diff. Eq. 71, 34 (1988); Ergod. Theor. Dynam. Sys. 8, 119 (1988). 
M. Grillakis, Comm. Pure Appl. Math. 41, 747 (1988); 43, 299 (1990). 

M. Grillakis, J. Shatah, and W. Strauss, J. Funct. Anal. 74, 160 (1987); ibid 94, 308 (1990). 

D. E. Pelinovsky, A.V. Buryak, and Yu.S. Kivshar, Phys. Rev. Lett. 75, 591 (1995); D.E. Pelinovsky, V.V. Afanasjev, and 
Yu.S. Kivshar, Phys. Rev. E 53, 1940 (1996). 



10 



[17] D.E. Pelinovsky and R.H.J. Grimshaw, In: Advances in Fluid Mechanics, Vol. 12, Eds. L. Debnath and S.R. Choudhury 

(Comp. Mech. Publ. Southampton, 1997), p. 246. 
[18] A.V. Buryak, Yu. S. Kivshar, and S. Trillo, Phys. Rev. Lett. 77, 5210 (1996); C. Etrich, U. Peschel, F. Lederer, and B.A. 

Malomed, Phys. Rev. E 55, 6155 (1997); A. De Rossi, C. Conti, and S. Trillo, Phys. Rev. Lett. 81, 85 (1998). 
[19] E.A. Ostrovskaya, Yu.S. Kivshar, D. Skryabin, and W. Firth, Phys. Rev. Lett. 83, 296 (1999). 
[20] D.V. Skryabin, Physica D 139, 186 (2000). 

[21] I.V.Barashenkov, D.E. Pelinovsky, and E.V. Zemlyanaya, Phys. Rev. Lett. 80, 5117 (1998). 

[22] Y.A. Li and K. Promislow, Physica D 124, 137 (1998); D. Mihalache, D. Mazilu, and L. Torner, Phys. Rev. Lett. 81, 4353 

(1998); M. Johansson and Yu.S. Kivshar, Phys. Rev. Lett. 82, 85 (1999). 
[23] L. Berge, Phys. Rev. E 58, 6606 (1998); O. Bang, L. Berge, and J.J. Rasmussen, Phys. Rev. E 59, 4600 (1999). 
[24] A.C. Yew, B. Sandstede, and C.K.R.T. Jones, Phys. Rev. E 61, 5886 (2000). 

[25] J. Yang and D.J. Benney, Stud. Appl. Math. 96, 111 (1996); J. Yang, Stud. Appl. Math. 98, 61 (1997); D.E. Pelinovsky 
and J. Yang, Stud. Appl. Math. (2000). 

Figure Captions 

Fig. 1: Eigenvalues 7 versus \i in the problem A(/j,)i/ = for N = 3: (a) a stable problem with no roots 

of 7(/i) for fi < 0, when p(U) = n(A) = 3; (b) an unstable problem with a single root of 7^) for fi < 0, when 
p(U) = 3 < n(A) = 4. 



11 




12 




13 



